Identification of molecular subtypes and prognostic signatures based on transient receptor potential channel-related genes to predict the prognostic risk of hepatocellular carcinoma: A review

Abnormal transient receptor potential (TRP) channel function interferes with intracellular calcium-based signaling and causes malignant phenotypes. However, the effects of TRP channel-related genes on hepatocellular carcinoma (HCC) remain unclear. This study aimed to identify HCC molecular subtypes and prognostic signatures based on TRP channel-related genes to predict prognostic risks. Unsupervised hierarchical clustering was applied to identify HCC molecular subtypes using the expression data of TRP channel-related genes. This was followed by a comparison of the clinical and immune microenvironment characteristics between the resulting subtypes. After screening for differentially expressed genes among subtypes, prognostic signatures were identified to construct risk score-based prognostic and nomogram models and predict HCC survival. Finally, tumor drug sensitivities were predicted and compared between the risk groups. Sixteen TRP channel-related genes that were differentially expressed between HCC and non-tumorous tissues were used to identify 2 subtypes. Cluster 1 had higher TRP scores, better survival status, and lower levels of clinical malignancy. Immune-related analyses also revealed higher infiltration of M1 macrophages and higher immune and stromal scores in Cluster 1 than in Cluster 2. After screening differentially expressed genes between subtypes, 6 prognostic signatures were identified to construct prognostic and nomogram models. The potential of these models to assess the prognostic risk of HCC was further validated. Furthermore, Cluster 1 was more distributed in the low-risk group, with higher drug sensitivities. Two HCC subtypes were identified, of which Cluster 1 was associated with a favorable prognosis. Prognostic signatures related to TRP channel genes and molecular subtypes can be used to predict HCC risk.


Introduction
Hepatocellular carcinoma (HCC) is a primary malignant hepatocyte tumor, accounting for approximately 80% of liver cancer cases worldwide. [1,2] According to annual projections by the World Health Organization, >1 million patients are expected to die of liver cancer by 2030. [3] The progression of cirrhosis is the biggest risk factor for HCC, and it is estimated that 90% of patients with underlying cirrhosis will eventually develop HCC. [4] Previous studies have demonstrated the beneficial chemopreventive effect of insulin-sensitizing drugs or statins against HCC occurrence. [5,6] Surgical resection, radiofrequency ablation, and liver transplantation are the first therapy choices for early HCC, with a 5-year survival rate of 70%. [7,8] During the last 15 years several aspects of HCC scenario have changed, as well as its management. [9,10] However, these therapies are ineffective in advanced stages (such as in metastatic HCC) where tumor resection is impossible, resulting in a 5-year survival rate as low as 2.5%. [11] A series of genetic events contribute to the pathogenesis and progression of HCC, 25% of which result from inducible mutations. These molecular alterations determine the abnormal proliferation of tumor cells and their subsequent invasion. [12] Therefore, biomarkers characterizing molecular genetic mechanisms are essential for assessing disease prognosis and predicting treatment responses.
Transient receptor potential (TRP) channels are a class of cation channels that perform signal transduction upon activation by changing the membrane potential or intracellular calcium (Ca2 + ) concentration [13] Ca 2+ signaling is crucial for regulating key cellular events, including gene transcription, movement and contraction, energy production, and channel control. [14] The mammalian TRP superfamily consists of 28 nonselective cation-permeable channels that can be divided into 6 subfamilies based on sequence homology. [15] In cancer, mutations in TRP channel-encoding genes generate TRP channels with abnormal functions that interfere with normal intracellular Ca 2+ distribution patterns, resulting in the dysregulation

Two molecular subtypes of hepatocellular carcinoma
were identified based on 16 transient receptor potential channel-related genes. 2. Cluster 1 was associated with a favorable prognosis and was more prevalent in the low-risk group. 3. Six prognostic signatures were used to calculate risk scores, which were defined as independent prognostic risk factors. 4. Risk score-based prognostic models can be used to stratify and predict the risk of HCC. of downstream effectors and enhanced cancer-specific pathological features. [16] For example, the expression of TRPV2 has been implicated in the drug-induced cytotoxicity and stemness of HCC. [17,18] In addition, the expression of TRPV2 increases as liver disease progresses from normal liver to chronic hepatitis and then to cirrhosis and may serve as a prognostic marker for patients with HCC. [19] Another TRP superfamily gene, TRPC7, is highly expressed in viral hepatitis B-associated HCC and may be a potential therapeutic target or diagnostic marker. [20] Therefore, genes encoding TRP channels may play crucial roles in the progression of HCC. However, the understanding of the role of all 28 channel-related genes is still incomplete, and systematic studies on their impact on the prognosis of HCC are lacking.
To this end, we analyzed 28 TRP channel genes from public expression data and established HCC molecular subtypes and prognostic signatures to stratify and predict the prognostic risk of HCC. We further constructed and validated a risk prognosis model and analyzed the resulting prognostic values reported for the immune microenvironment characteristics and drug sensitivity of TRP channel genes to comprehensively elucidate the potential role of these genes in HCC.

Data sourcing and preprocessing
Standardized gene expression profiles and clinical follow-up data of patients with HCC were obtained from The Cancer Genome Atlas (TCGA) database. After removing samples with no overall survival data, 365 HCC samples with prognostic information were retained. The microarray dataset GSE14520 [21] downloaded from the Gene Expression Omnibus database was defined as an independent external validation cohort. GSE14520 is present on the sequencing platform of the GPL3921 [HT_HG-U133A] Affymetrix HT Human Genome U133A Array and contains 221 HCC samples with valid prognostic data.

Subtype identification based on TRP channel genes
Data on the 28 TRP channel genes used in this study are provided in a relevant published article. [22] We compared expression levels between HCC and non-tumorous samples and used Pearson correlation analysis to assess potential associations. We identified the corresponding HCC subtypes based on TRP channel-related genes with significant differences in expression, which were subjected to unsupervised hierarchical clustering using R 3.6.1 ConsensusClusterPlus Version 1.54.0 (http://www.bioconductor.org/packages/release/bioc/ html/ConsensusClusterPlus.html). The optimal k value ranged from 2 to 6. The TRP score in each sample was determined by enrichment scores computed by the gene set variation analysis algorithm using the R gene set variation analysis package version 1.36.2 (http://bioconductor.org/packages/release/bioc/ html/GSVA.html) [23] and subsequently compared among subtypes using the Wilcoxon test to verify the rationality of HCC subtypes.

Association analysis between HCC subtypes and clinical features
To evaluate the correlation between prognostic survival and different HCC subtypes, Kaplan-Meier (KM) curves were generated using the R survival package Version 2.41-1 (http:// bioconductor.org/packages/survivalr/). [24] The clinical characteristics of patients with HCC (including age, sex, tumor grade, pathological stage, pathological TNM classification, and exposure to radiation therapy) were analyzed. Furthermore, the relationships between subtypes and these clinical features were assessed using the chi-square test, and statistical significance was set at P < .05.

Comparison of immune microenvironment among subtypes
In this study, we applied 2 algorithms to evaluate the immune microenvironment status of HCC samples and compared immune infiltration among different HCC subtypes using the Wilcoxon signed-rank test. For immune infiltration analysis, CIBERSORT (httRiskscore://cibersort. stanford.edu/index.php) [25] was used to compute the proportions of 22 types of immune cells, whereas ESTIMATE [26] was used to estimate the stromal and immune scores of the tumor samples.

Screening of differentially expressed genes (DEGs) among HCC subtypes
To observe the possible existence of different molecular mechanisms underlying HCC subtypes, DEGs between these subtypes were screened using linear regression and empirical Bayesian methods provided by the limma package, version 3.10.3 (http:// www.bioconductor.org/packages/2.9/bioc/html/limma.html). [27] The obtained P values were adjusted for multiple tests using the Benjamini-Hochberg method, and genes characterized by an adjusted P value < 0.05 and |log 2 fold change| > 1 were considered to have significant differential expression.

Identification of prognostic signatures
Based on the obtained DEGs, univariate Cox regression analysis was used to screen the expression levels of genes significantly associated with survival with a set threshold of P < .01. The least absolute shrinkage and selection operator algorithm was then used to define an optimal lambda value, followed by subsequent assessment of prognostic signatures (10-fold cross-validation) using the R package for the LARS algorithm version 1.2 (https:// cran.r-project.org/web/packages/lars/index.html).

Generation and validation of the risk score-based prognostic model
Furthermore, the prognostic gene signatures were subjected to stepwise Cox regression analysis using the R survminer package version 0.4.9 (https://cran.rstudio.com/web/packages/ survminer/index.html) to establish a risk score-based prognostic model. The risk score was calculated using the following formula: where β is the regression coefficient, h 0 (t) is the baseline risk function, and h(t, X) is the risk function associated with X Figure 3. Comparison of immune cell infiltration between cluster 1 and cluster 2 using CIBERSORT and ESTIMATE algorithms. (A) Infiltration differences in 22 types of immune cells between the 2 clusters were evaluated using the CIBERSORT algorithm. *P < .05, **P < .01, ***P < .001. (B) The violin plot depicts the difference in stromal and immune scores between clusters 1 and 2, which were calculated using the ESTIMATE algorithm. www.md-journal.com (covariate) at time t. We further calculated the risk scores of each sample in TCGA and GSE14520 cohorts and grouped the samples according to their median values. The samples were divided into high-risk prognosis (HRP) and low-risk prognosis (LRP) groups, and a KM curve was created to estimate the difference in the actual prognosis.

Analysis of prognostic independence to construct a nomogram model
We assessed clinical characteristics using univariate and multivariate Cox regression analyses to screen for independent prognostic factors at a threshold value of P < .05. The R.rms package version 5.1-2 (https://cran.r-project.org/web/packages/rms/index.html) [28] was used to construct a nomogram model to predict the survival probability of patients with HCC.

Drug sensitivity analysis
The Genomics of Drug Sensitivity in Cancer database was used to estimate the sensitivity of each sample to chemotherapy drugs, followed by quantification of 50% inhibitory concentration (IC50) values using the R pRRophetic package (https://github.com/paulgeeleher/pRRophetic). [29] The IC50 values of 138 chemotherapy drugs were compared between the HRP and LRP groups using the Wilcoxon signed-rank test.

Association analysis of prognostic risk and HCC subtype
To observe the relationship between HCC subtypes and risk groups, we compared the distribution differences of subtypes between the HRP and LRP groups. The results were visualized using a bar chart obtained with the R. galluvial package version 0.12.3 (https://CRAN.R-project.org/package=ggalluvial). [30] 3. Results

Two HCC subtypes were identified based on the analysis of TRP channel-related genes
To analyze the expression patterns of the 28 TRP channel-related genes, their expression profiles were compared between HCC and non-tumorous samples. Sixteen TRP channel-related genes showed significant differences in expression levels between the 2 groups (Fig. 1A). Correlations between the expression of these genes were also analyzed; the results suggested that TRPV4 and TRPV6 had the strongest positive correlation, whereas TRPM2 and TRPM7 had the strongest negative correlation (Fig. 1B). The 16 key TRP-related genes were then incorporated into the unsupervised cluster analysis, and the optimal subtype number was defined as k = 2 (Fig. 1C). The proportion of ambiguous clustering was then determined to verify the stability of the clustering results, and the optimal value of k = 2 was confirmed (Fig. 1D). The TRP score of each sample was calculated and compared between the 2 clusters. The TRP score of Cluster 1 was significantly higher than that of Cluster 2 (Fig. 1E). Survival analysis also indicated a poor prognosis in patients with a low TRP score (Fig. 1F). These results highlight that a higher TRP value, as exhibited by Cluster 1, may be a favorable prognostic factor.

Differences in clinical features and prognosis between subtypes
To further explore the relationship between clinical characteristics and HCC subtypes, we conducted a survival analysis to compare the survival differences between Clusters 1 and 2 ( Fig. 2A). The KM curve suggested a better prognosis for patients in Cluster 1, which is consistent with previous results. In addition, the distribution of expression of the 16 key TRP-related genes associated with different clinical features was visualized using a heatmap (Fig. 2B). Comparison among subtypes revealed that the distribution of Cluster 2 was significantly associated with the later period of pathologic T, pathologic N, and pathologic stage (Fig. 2C). This suggests that Cluster 2 might possess a prognostic risk value, indicating the deterioration of clinical manifestations.

Analysis of immune microenvironment characteristics based on HCC subtype
Based on the expression matrix of HCC samples, the CIBERSORT algorithm was employed to calculate the fractions of immune cells in each sample. Using a threshold value of P < .05, 10 types of immune cells were found to infiltrate significantly differently between Clusters 1 and 2 (Fig. 3A). Among these, M1 and M2 macrophages exhibited the most significant difference between the 2 groups; their infiltration levels were higher in Cluster 1 than in Cluster 2. The ESTIMATE algorithm was used to compute the stromal and immune scores for each sample. Compared to Cluster 2, Cluster 1 presented higher immune and stromal scores (Fig. 3B), which also indicated a higher level of non-tumor cell infiltration.

Optimal prognostic signatures were screened based on DEGs between HCC subtypes
To further explore the differences in molecular regulatory mechanisms among the subtypes, we performed a differential analysis of gene expression levels between Clusters 1 and 2. The patterns of the DEGs between the 2 groups are shown in Figure 4A. A total of 148 upregulated and 530 downregulated DEGs were identified according to the set thresholds (Fig. 4B). Univariate Cox regression analysis was performed to select the 212 DEGs that were significantly associated with prognosis. The least absolute shrinkage and selection operator algorithm was used to screen for optimal candidates, and 16 key genes were selected at lambda = 16 (Fig. 4C). Finally, stepwise Cox regression was used to define the optimal gene set, which constituted the basis  SPP1, BRSK2, OGN, PPARGC1A, and FTCD) (Fig. 4D).

Prognostic model construction and efficacy verification
Based on the regression coefficients of the 6 prognostic signatures and their levels in the TCGA dataset, a risk score-based prognostic model was constructed. The distributions of the risk score, survival time, and 6 prognostic signatures for each sample are shown in Figure 5A. These results suggest that the risk of death increases with the risk score. Furthermore, the expression of MMP1, SPP1, and BRSK2 was higher in the HRP group than in the LRP group. The KM curve confirmed that the actual survival of patients in the LRP group was significantly higher than that of patients in the HRP group (Fig. 5B). A time-dependent receiver operating characteristic curve was created to assess the specificity and sensitivity of the model in predicting prognostic risks; the area under the curve (AUC) of the 1-, 3-, and 5-year prediction were 0.749, 0.742, and 0.712, respectively (Fig. 5C), which represented a high predictive performance of the model. The performance of the model based on the 6 prognostic signatures was analyzed using the GSE14520 validation dataset. The results suggested that with an increase in the risk score, patient survival time and gene expression levels tended to change (Fig. 5D). The patients in the HRP group also had a significantly lower likelihood of survival (Fig. 5E). In addition, all AUC values were >0.6, indicating that the prognostic model based on the validation cohort had the potential to predict prognostic risk (Fig. 5F).

Screening of independent prognostic factors to construct a nomogram survival prediction model
We applied univariate Cox regression analysis for the clinical characteristics of all HCC samples in the TCGA dataset and selected factors with P < .05 to further conduct the multivariate Cox regression analysis ( Fig. 6A and B). Then, the factors "pathological T" and "risk groups" were identified with prognostic independence to establish a nomogram model (Fig. 6C). The calibration curves confirmed that the predicted overall survival at 1, 3, and 5 years using the nomogram model fit well with the actual values (Fig. 6D). In the KM curve, patients in the HRP group (as predicted by the nomogram model) also showed a significantly unfavorable prognosis compared to those in the LRP group (Fig. 6E). Additionally, the AUC of the 1-, 3-, and 5-year receiver operating characteristic curves were 0.761, 0.739, and 0.703, respectively, confirming the predictive reliability of the nomogram model (Fig. 6F).

Associations between prognostic risk with HCC subtypes and drug sensitivity
To observe the relationship between the HCC subtypes and risk groups, we compared the distribution of the 2 clusters between the HRP and LRP groups. The IC50 values of all drugs were significantly higher in the LRP group than in the HRP group (Fig. 7A), indicating that the high-risk group has higher drug sensitivities. Finally, 64 of the 138 chemotherapeutic drugs showed significantly different IC50 values between the HRP and LRP groups. Using box plots (Fig. 7B), we listed the differences in IC50 values between the 2 groups for the 6 chemotherapy drugs frequently used for HCC. The results showed that the IC50 values of these drugs were lower in the high-risk group than in the lower-risk group, indicating that differences in drug sensitivity may be a source of prognostic risk for patients with HCC.

Discussion
In this study, 16 of the 28 TRP channel-related genes were differentially expressed between HCC and non-tumorous sample tissues. Based on these key TRP channel-related genes, we identified 2 molecular HCC subtypes, of which Cluster 1 had a higher TRP score and possibly a more favorable prognosis. Analysis of prognostic factors and clinical features confirmed the prognostic characteristics of Cluster 1. Compared with patients in Cluster 2, more patients in Cluster 1 were distributed in the early stage of HCC and were characterized by a low malignant level of the disease. After screening the DEGs between Clusters 1 and 2 and constructing the prognostic model, we found that prognostic signatures based on the HCC subtype could predict prognostic risks effectively. Moreover, patients in Cluster 1 were more likely to be distributed in the low-risk group than those in Cluster 2. These results corroborate each other, demonstrating the reliability of subtype-based clustering and the successful construction of the prognostic model and confirming the important role of the 16 TRP channel-related genes in predicting the prognosis of HCC. Taken together, the prognostic signatures related to TRP channel genes and molecular subtypes obtained in this study can be used in clinical applications to predict HCC risk.
Identification of the immune microenvironment characteristics of the 2 subtypes revealed that the most significant difference between Clusters 1 and 2 was in the infiltration levels of M1 macrophages. Macrophages play a vital role in many pathophysiological processes such as inflammation, tissue repair, and metabolism. The polarization of M1/M2 macrophages may affect immune escape, carcinogenesis, and metastasis in HCC [31,32] However, in our study, no difference between clusters was detected in terms of the infiltration level of M2 macrophages. However, a relevant study identified DEGs between high-and low-infiltrating M1 macrophages as predictors of HCC prognosis [33] The activation of Notch signaling in HCC can mediate the differentiation of macrophages into M1 macrophages, thereby promoting inflammation and anti-tumor activity. [34] In this study, a significantly increased infiltration of M1 macrophages was observed in Cluster 1, indicating the activation of the inflammatory response and anti-tumor activity against HCC. This is why patients in Cluster 1 showed better clinical outcomes. M1 macrophages play an anti-tumor role by inhibiting tumor cell proliferation and metastasis, [35] which is consistent with our findings. Another finding is that the activation of the TRP family-related gene, TRPC5 (which in this study was highly expressed in paracarcinoma tissues), may inhibit macrophage differentiation by regulating the Akt/IκB/NF-κB signaling pathway. [36] Therefore, we further speculated that TRPC5 could cause increased infiltration of M1 macrophages, thereby reducing the invasiveness of tumor cells, inhibiting their malignant phenotype, and promoting a favorable prognosis for patients with HCC. However, the regulatory mechanism of M1 macrophages and the role of M1/M2 polarization in HCC remain to be experimentally explored.
Considering the differences in overall survival and clinical phenotypes between the 2 HCC subtypes, we further analyzed the DEGs to explore the underlying molecular regulatory differences between Clusters 1 and 2 and screened prognostic signatures for the construction of a prognostic model and prognostic risk prediction. Among these, FTCD, PPARGC1A, OGN, BRSK2, SPP1, and MMP1 have important prognostic values and can be used to define prognostic risks. Specifically, FTCD, PPARGC1A, and OGN are protective factors against HCC and are associated with a better prognosis. Chen et al confirmed our findings and proposed that low expression of FTCD is associated with poor prognosis and an aggressive tumor phenotype in HCC. [37] The sensitivity and specificity of FTCD expression levels in differentiating normal or liver cirrhosis from early well-differentiated HCCs are 90% and 86.7%, respectively. [38] Conversely, BRSK2, SPP1, and MMP1 expression levels were associated with a higher risk of poor prognosis in this study. Li et al further demonstrated that BRSK2 is positively associated with cancer status, has a prognostic risk for HCC, and may be involved in the regulatory mechanism of m6A methylation. [39] SPP1 is also a signature gene in HCC tissues that enhances the proliferation of tumor cells and is closely related to tumor progression. [40] More importantly, SPP1 mediates the crosstalk between HCC cells and macrophages and triggers the polarization of macrophages towards the M2 phenotype. [41] However, the association between these genes and the TRP family has rarely been studied, except for MMP1. MMP1 is upregulated in patients with HCC and is associated with a poor prognosis. [42] Pharmacological inhibition and gene silencing of TRPA1 downregulates MMP1 production in osteoarthritis, [43] whereas activation of TRPV1 increases the expression of MMP1. [44] However, in the present study, specific levels of TRPA1 and TRPV1 were not detected in HCC tissues. Therefore, the relationship between the expression of these prognostic gene signatures and the TRP channel genes requires further confirmation.
There were several limitations in this study. First, there was no analysis of the potential associations between key TRP genes, M1 macrophage infiltration, and prognostic signatures. Second, no experimental verification of the diagnostic and prognostic values of key genes was performed. Third, all analysis was based on public databases, and no functional assay was conducted. Hence, in subsequent studies, we will focus on the interaction between M1/M2 macrophage polarization and TRP channel signaling, as well as on their actual impact on the prognosis of HCC.

Conclusion
In this study, 2 HCC subtypes were identified based on key TRP channel-related genes. Patients in Cluster 1 were prone to a favorable prognosis and a lower level of tumor malignancy. Based on differential gene screening between Clusters 1 and 2, and prognostic signature identification, we constructed prognostic and nomogram models to predict and stratify the prognostic risks of HCC. The validation results confirm the robust prediction performance of the models. However, subsequent experimental verification is required.